Spatial Statistics Research - Predicting crimes in Los Angeles

Number of Crimes VS Average Household Size

The spatial data delivers more detailed information on certain underlying relationships. Under the hypothesis stated in the introduction, one would expect that the more a circle tends to red (i.e. the higher the average household size), the more an area tends to blue (the higher the number of crimes). This is however not the case, except for three areas north-west of Santa Monica.

## Mapping Number of Crimes VS Average Household Size

centroids <- as.data.frame(gCentroid(zCrimes,byid=TRUE))

qpal <- colorQuantile("YlGnBu", zCrimes$crimes, n = 5)
qHH <- colorQuantile("YlOrRd", zCrimes$Average.Household.Size, n = 5)

leaflet(zCrimes) %>% addPolygons(weight = 1, smoothFactor = 0.5,
    opacity = 1.0, fillOpacity = 0.5,
    color = ~qpal(crimes),
    highlightOptions = highlightOptions(color = "white", weight = 2,
      bringToFront = TRUE)) %>% 
                addCircles(lng = ~centroids$x, lat = ~centroids$y, weight = 1, color = ~qHH(Average.Household.Size),
                  radius = 800, popup = ~zCrimes$zip, opacity = 0.9, fillOpacity = 0.8 ) %>% 
                addTiles() %>%
    addLegend(pal = qpal, values = ~crimes, opacity = 1) %>%
    addLegend(pal = qHH, values = ~Average.Household.Size, opacity = 1)

Number of Crimes VS Median Age

Similarily, for median age against crime rates no particular relationship can be read from the plot.

## Mapping Number of Crimes VS Median Age

qMA <- colorQuantile("YlOrRd", zCrimes$Median.Age, n = 5)

leaflet(zCrimes) %>% addPolygons(weight = 1, smoothFactor = 0.5,
    opacity = 1.0, fillOpacity = 0.5,
    color = ~qpal(crimes),
    highlightOptions = highlightOptions(color = "white", weight = 2,
      bringToFront = TRUE)) %>% 
                addCircles(lng = ~centroids$x, lat = ~centroids$y, weight = 1, color = ~qMA(Median.Age),
                  radius = 800, popup = ~zCrimes$zip, opacity = 0.9, fillOpacity = 0.8 ) %>% 
                addTiles() %>%
    addLegend(pal = qpal, values = ~crimes, opacity = 1) %>%
    addLegend(pal = qMA, values = ~Median.Age, opacity = 1)

Ordinary Linear Regression

Below is a regression analysis for the two independent variables (average household size and median age). It is deceiving in terms of explained variance and of the explanatory value of the independent variable median age. Average household size is significantly different from zero at a 5% confidence level though: a bigger household would imply more crimes as stated in the hypothesis in the introduction. Generally a high R-squared is a sign of possible multicollinearity. Although R-squared is particularily low here, one might still need to perform formal tests in order to assess the absence of multicollinearity.

regHH <- lm(log(crimes) ~ Average.Household.Size + Median.Age, data = zCrimes)
stargazer(regHH, no.space = T, type = "html")

# Multicollinearity

mc <- omcdiag(zCrimes@data[colnames(zCrimes@data) %in% c("Average.Household.Size", "Median.Age")], log(zCrimes$crimes))

kable(mc$odiags, "html", booktabs = T, caption = "multicollinearity tests") %>%
  kable_styling("striped")

# HETEROSKEDASTICITY

# ols_test_bartlett(as.data.frame(zCrimes@data$"Average.Household.Size"), as.data.frame(zCrimes@data$"Median.Age"))
bp <- ols_test_breusch_pagan(regHH)

m <- matrix(c(bp$bp, bp$p, 1))
rownames(m) <- c("Chi2", "Prob > Chi2", "DF")
kable(m, "html", booktabs = T, row.names = T, caption = "Breusch-Pagan test") %>% kable_styling("striped")

# Spatial Dependence

spatmatrix <- poly2nb(zCrimes)
w <- nb2listw(spatmatrix, zero.policy = T)

# Moran's I
moran <- lm.morantest(regHH, listw = w, zero.policy = T)# Moran’s I test for residual spatial autocorrelation
moranStat <- matrix(c(moran$statistic, moran$p.value, 1))
rownames(moranStat) <- c("statistic", "p.value", "parameter")
colnames(moranStat) <- "Moran's I (error)"

# Lagrange Multiplier test
LM <- lm.LMtests(regHH, listw = w, zero.policy = T, test = c("LMerr", "RLMerr", "LMlag", "RLMlag", "SARMA")) 

LMrows <- c("statistic", "p.value", "parameter")
m <- matrix(NA, nrow = 3, ncol = length(LM))
rownames(m) <- LMrows
colnames(m) <- c("Lagrange Multiplier (error)", "Robust LM (error)", "Lagrange Multiplier (lag)", "Robust LM (lag)", "Lagrange Multiplier (SARMA)")
# colnames(m) <- 
  
# everything in one table
for(i in 1:length(LM)){
  r <- 1
  for(j in LMrows){
    m[r, i] <- LM[[i]][[j]]
    r <- r + 1
  }
}
m <- cbind(moranStat, m)
kable(m, "html", booktabs = T, row.names = T, caption = "spatial dependence tests")  %>%
  kable_styling("striped")

# Residuals

#extract residuals
zCrimes@data$lmRes[!is.na(zCrimes$crimes)] <- resid(regHH) #residual lm

#view residuals for linear model
qpal <- colorQuantile("OrRd", zCrimes@data$lmRes, n=9) 

leaflet(zCrimes) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(lmRes)
  ) %>%
  addTiles()

OLS with lag or error term to control for spatial autocorrelation

In order to deal with spatial autocorrelation, an OLS with spatial lag and an OLS with spatial errors are fitted. Model selection will be based on AIC.

#run a spatial lag model
lagModel <- lagsarlm(log(crimes) ~ Average.Household.Size + Median.Age, data = zCrimes, w,
                  zero.policy = T)
stargazer(lagModel, no.space = T, type = "html")

#note that rho is listed below other model coefs.

#run a spatial error model
errModel <- errorsarlm(log(crimes) ~ Average.Household.Size + Median.Age, data = zCrimes, w, zero.policy = T)
stargazer(errModel, no.space = T, type = "html")

#extract residuals for lag model
zCrimes@data$lagRes[!is.na(zCrimes$crimes)] <- resid(lagModel) #residuals lag

#view residuals for lag model
qpal<-colorQuantile("OrRd", zCrimes@data$lagRes, n=9) 

title <- tags$div(tag.map.title, HTML("spatial lag model residuals"))  

leaflet(zCrimes) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(lagRes)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")

#extract residuals for err model
zCrimes@data$errRes[!is.na(zCrimes$crimes)] <- resid(errModel) #residual err

#compare with residuals for err model
qpal<-colorQuantile("OrRd", zCrimes@data$errRes, n=9) 

title <- tags$div(tag.map.title, HTML("spatial error model residuals"))  

leaflet(zCrimes) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(errRes)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")

Geographically Weighted Regression Model

Another way of dealing with the intrinsic spatial interactions is to fit a geographically weighted regression model.

centr <- gCentroid(zCrimes, byid = TRUE)  #byid tells function to add centroids per block group

# create SpatialPointsDataFrame to add all data from chi.poly@data to centr@data
centr <- SpatialPointsDataFrame(centr, data= zCrimes@data) 


#Begin GWR prep: First, calculate bandwidth

#calculate optimal bandwidth for model (uses cross validation to find kernel bandwidth that
#generates best model...e.g.-kernel that minimizes RSS)

GWRbandwidth <- gwr.sel(log(crimes) ~ Average.Household.Size + Median.Age, data=centr, verbose = F, 
                        coords=coordinates(centr),adapt=T)# coordinates() extracts x, y coordinates
#adapt tells function to calculate adaptive bandwidth

#run the gwr model
gwr.model <- gwr(log(crimes) ~ Average.Household.Size + Median.Age, 
                 data=subset(centr, !is.na(crimes)), coords=coordinates, adapt=GWRbandwidth,
                 hatmatrix=TRUE, se.fit=TRUE)
#print the results of the model

#extract results for each variable
results<-as.data.frame(gwr.model$SDF)

#attach coefficients to original dataframe
centr$Median.AgeCoef[!is.na(zCrimes$crimes)]  <- results$Median.Age
centr$Average.Household.SizeCoef[!is.na(zCrimes$crimes)]  <- results$Average.Household.Size

#now plot the various GWR coefficients

#Does the median age influence violent crime uniformly across LA?

title <- tags$div(tag.map.title, HTML("median age influence on violent crime"))  

qpal<-colorQuantile("OrRd", centr@data$Median.Age, n=9) 

leaflet(centr) %>%
  addCircleMarkers(radius=3,color = ~qpal(centr@data$Median.Age)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")

The plot above shows where median age influence violent crime across LA. It seems that patches of strongly influenced neighborhoods exist, implying that the influence is not uniform accross the city.

#Does the average household size influence violent crime uniformly across LA?
title <- tags$div(tag.map.title, HTML("average household size influence on violent crime"))  

qpal<-colorQuantile("OrRd", centr@data$Average.Household.Size, n=9) 

leaflet(centr) %>%
  addCircleMarkers(radius=3,color = ~qpal(centr@data$Average.Household.Size)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")

The plot above shows where average household size influence violent crime across LA. Similarily, it seems that patches of strongly influenced neighborhoods exist. A clearer pattern is distinguishable: the outskirts of the city are more influenced by average household size. This might be due to the fact that poorer people live in the outskirts and that a bigger household has a more significant impact on the ability to foster and educate children.